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Abstract. We investigate both numerically and analytically the effect of strong disorder on the large scale 
properties of the hyperbolic equations for stresses proposed in Jl],[|. The physical mechanism that we 
model is the local splitting of the force chains (the characteristics of the hyperbolic equation) by packing 
defects. In analogy with the theory of light diffusion in a turbid medium, we propose a Boltzmann-like 
equation to describe these processes. We show that, for isotropic packings, the resulting large scale effective 
equations for the stresses have exactly the same structure as those of an elastic body, despite the fact that 
no displacement field needs to be introduced at all. Correspondingly, the response function evolves from a 
two peak structure at short scales to a broad hump at large scales. We find, however, that the Poisson ratio 
is anomalously large and incompatible with classical elasticity theory that requires the reference state to 
be thermodynamically stable. 



PACS. 05.40.-a Fluctuation phenomena, random processes, noise, and Brownian motion 
sandpiles; granular compaction - 83.70.Fn Granular solids 



45.70.Cc Static 



1 Introduction 

The mechanics of assemblies of hard, cohesionless grains 
has been the subject of rather controversial studies in the 
past few years ]|,[|||||. A well developed approach, in 
particular for engineering applications, is based on the as- 
sumption that granular materials behave, on large length 
scales, as elastic (or elastoplastic) bodies. The microscopic 
justification of this assumption for undeformable cohesion- 
less particles is however far from obvious. The conven- 
tional theory of elastic bodies starts from the identification 
of a reference state from which deformations are defined; 
the free-energy of the deformed body is then expanded in 
powers of the strain tensor . The definition of a smooth 
deformation field and of its energy are both problematic 
in the context of hard grains || . Recent numerical simu- 
lations furthermore suggest that a limiting curve relating 
stresses and deformations for large enough volumes might 
not even exist || . 

The absence of any obvious deformation field from 
which the stress tensor may be constructed has motivated 
an alternative, 'stress-only' approach Jl],||,|]]|] . The basic 
tenet of these theories is that in equilibrium, some (his- 
tory dependent) relations between the components of the 
stress tensor are established. A well known relation of this 
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type arises from the assumption that the material is ev- 
erywhere on the verge of plastic failure, leading to a Mohr- 
Coulomb relation between the stress components pLOj ] . A 
simpler relation, based both on symmetry arguments and 
on the consideration of simple rules for the transfer of 
stresses between adjacent grains, is a local proportional- 
ity between the diagonal components of the stress tensor 
[01. For example, in two dimensions, this relation reads 
a zz — Cq<j xx , where Co is a constant. This equation is a 
local version of the hypothesis made by Janssen in his 
famous theory for stresses in silos The conse- 

quence of this apparently innocuous assumption is that 
stresses obey an hyperbolic equation, as compared to the 
elliptic equations encountered in elasticity theory. This 
means that stresses 'propagate' along lines: the charac- 
teristics of this hyperbolic equation were argued to be the 
mathematical transcription of the force chains that are 
well known to exist in granular materials |l2| . More elab- 
orate versions of the above closure relation between the 
components of the stress tensor have been investigated in 
|U^], and have been shown to reproduce both th e p res- 
sure 'dip' underneath the apex of a sandpilc [|l3Ul4|| , or 
the pressure profile inside silos JllJ . 

A rather clear-cut difference between hyperbolic and 
elliptic equations lies in the structure of the propagator, 
i.e. the response of the stress to a localized force applied 
at the top of a layer of granular matter of height h. For 
an hyperbolic equation, the pressure profile at the bottom 
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of the layer is localized around two peaks (in two dimen- 
sions) or a circular annulus (in three dimensions), centred 
underneath the point where the force is applied For 
an elliptic equation, by contrast, the response function is a 
belly shaped curve of width oc h that reaches its maximum 
precisely underneath the point where the force is applied. 
A similar single peak structure, but of width oc Vh, is also 
predicted by scalar models |lq| , where the stress obeys a 
(parabolic) diffusion equation. 

This response function has only very recently been 
carefully measured experimentally. For well ordered pack- 
ings, or for frictionless beads, the two peak structure is 
rather convincingly observed both in two and three di- 
mensio ns |p^ ,p7|, and has been further justified theoreti- 
cally in [|18|. However, as soon as the packing is disordered 
(for example by considering mixtures of grains of different 
sizes), the two peak structure disappears on large length 
scales and is replaced by a one peak, elastic like, response 
function p9| , ^6p0| ], with a width scaling as h. This single 
hump response function has been also observed in numer- 
ical simulations of disordered packing in ^T|. Finally, a 
purely diffusive response function (scaling as yfh) was re- 
ported in |22]] , but on a very special 'brick' packing. 

The aim of this paper is to investigate both numeri- 
cally and theoretically the consequences of strong disorder 
on the large scale properties of the locally hyperbolic equa- 
tions proposed in fjjffl. We propose and simulate numer- 
ically a simple model where force chains (characteristics) 
are 'scattered' in a random way by defects. The response 
function is seen to evolve, as a function of the height h of 
the layer, from a two peak structure at small scales to a 
one peak elastic like response at large scales. We then pro- 
pose a Boltzmann-like formalism to describe the statistics 
of this force chain scattering process. We show that the 
large scale effective equations for the stresses become ellip- 
tic, and, rather surprisingly, have exactly the same struc- 
ture as the equations governing the stresses in an elastic 
body. This is quite striking because no displacement field 
needs to be introduced at all. (This is actually also true 
of our numerical model of where only forces are consid- 
ered) . Interestingly, it is the very existence of force chains 
that allows one to define a local vectorial 'order parame- 
ter' (that we call J below) which formally plays the role 
of the displacement in elasticity theory. The importance 
of force chains for the mechanics of granular materials 
was emphasized in p3| , [j"2"[ ; the present work suggests that 
their presence might be(somewhat paradoxically) crucial 
to allow for a large scale elastic like behaviour. 

The Poisson ratio v that appears in our elastic-like 
equations depends solely on the statistical properties of 
the force chain splitting process. Interestingly, this Pois- 
son ratio violates the bound v < 1/2 (in three dimensions) 
that comes from the fact that the elastic energy of defor- 
mations is a positive quadratic form of the strain tensor. 
This clearly shows that although the final equations have 
the same formal structure as those of elasticity theory, 
granular materials do not behave as standard elastic bod- 
ies, even for statistically isotropic packings. 




Fig. 1. Picture of the force chains in a two-dimensional system 
of grains subject to a vertical force imposed at the middle of 
the top surface |^(|. Such a picture has been obtained with 
birefringent grains between perpendicular polarizers. 



2 The force splitting model: definition and 
numerical simulations 

Pictures of photoelastic grains under compression clearly 
reveal the existence of linear force chains which tend to 
split upon meeting vacancies or packing defects: see figure 
|l| |2^,^5|. This motivates the following model to compute 
the response function G(x, h) to a force imposed at the top 
of the layer (x = 0, z = 0), and measured at x and z = h. 
For a perfectly ordered packing, we want this response to 
be made of two sharp peaks centered at x — icoft., where 
Co = tan$o is the tangent of the angle of the direction 
along which forces are transmitted in regular arrangement 
of grains. For a compact triangular packing, this angle is 
6*o = 30°, which is the value that we choose in the following 
simulations. Therefore we start with two force chains at 
angle ±9q to the vertical - the initial point x — z = 
itself can been seen as a 'defect', see below. 

For disordered packing, we randomly distribute non 
overlapping defects of radius a — 1 (this fixes the unit 
length) with a density p c i- If one of the force chains meets 
a defect, we split it into two new ones at random angle, 
which then propagate until another defect (or the bound- 
ary) is reached, an so on. 

More precisely, a chain carrying a force / in the direc- 
tion n splits when meeting a defect into two forces /i, 
in the directions ni, ri2. The two angles 9\ and 9i (be- 
tween n and ni and 112 respectively) are uniformly chosen 
between and some maximum splitting angle 9 m ■ The lo- 
cal mechanical equilibrium imposes that the intensities f± 
and /2 are such that /n = /ini + / 2 n 2 - see figure [2] (a). 
Sometimes, two (or more) force chains meet at the same 
defect. In this case, we make them all merge together - 
see figure |^ (b) . It is important to note that the positions 
of the defects are fixed before starting the computation of 
the forces. This idea of a frozen disorder is consistent with 
the experimental observation that when the local overload 
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Fig. 2. The disorder of the packing is modelled by the presence 
of 'defects' where force chains can split at random angles (a) 
('A processes') - or merge (b) ('Y processes'). The force chains 
have no width, but the defects have a finite radius a — 1. At 
each defect, forces are balanced: /n = /ini + /2ri2 . 



is added on the top of the system, the forces are transmit- 
ted along the chains originally created during the building 
of the packing. In other words, as long as the applied force 
is not too large and compatible with the pre-existing net- 
work of force chains, the geometry of the packing, and in 
particular the contacts between grains, remains the same. 

The network of forces is built in a hierarchical man- 
ner: each defect involved in this network has one or more 
'parents' and one or more 'children'. The most common 
situation is when a defect has one parent and two children. 
We qualify such defects to be of 'A type'. Defects of 'type 
Y' have two (or very rarely three or more) parents but 
one child only. Finally, some defects are not involved in 
the network of forces. The type of the defects is not fixed 
a priori. At the beginning of the calculation, all defects are 
'free' and become progressively of A type when reached by 
a force chain. When a force chain reaches a defect which is 
already of A type (or less likely of Y type), all force chains 
originally descending from this defect are removed, and it 
becomes of Y type. Note that this procedure does not al- 
low cycles, i.e. defects cannot be their own ancestors - or 
descendants [] We checked that the final network of force 
chains is independent of the order in which defects are 
chosen at each step. Eventually, all the force chains reach 
a boundary of the system, and the calculation stops. All 
four walls are perfect absorbing boundaries. 

A further technical difficulty is the presence of very 
small forces. The deeper the layer, the larger the num- 
ber of force chains that carry extremely weak forces. In 
order to keep the computation time to reasonable values, 
the total number of force chains cannot grow indefinitely. 
We therefore introduced a lower cut-off, below which force 
chains are no longer taken into account. We checked that 
for small enough cut-off our results are independent of this 
cut-off. 



Even if cycles were to be allowed, the problem would con- 
verge very quickly - there is no sense in which allowing cycles 
would modify, except infinitesimally, the results obtained. 
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Fig. 3. Snapshot of the network of force chains obtained within 
the AY model. The dots represent the defects of the granular 
packing. The lines are bolder when the forces are larger - very 
small forces (less than 0.01) have not been plotted. In this 
picture we set the height to h — 200, the width to 2w = 300, 
the density of defects to pd = 0.028 and the maximum splitting 
angle to 8 m = 30°. 



With these rules, realistic force networks can be cre- 
ated - see figure |. After averaging over many statistically 
identical samples, one can obtain stress profiles for differ- 
ent heights h. Figure shows that, as h increases, the ver- 
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Fig. 4. Averaged vertical stress response function for differ- 
ent values of h. At small scale, the response function has a 
double peak shape which is characteristic of a locally hyper- 
bolic behaviour. For larger heights, these two peaks merge into 
a single broad peak, comparable to an elastic-like response. 
These curves have been obtained with 2w = 500, pd = 0.028, 
6m = 30°, and averaged over 1000 samples. 



tical pressure response profile evolves continuously from 
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Fig. 5. For large heights, the width of the response function 
grows like h. The inset shows that its shape is not gaussian and 
actually rather close to a purely elastic response, shown as the 
plain line. These curves have been obtained with 2w — 2400, 
pd = 0.028, 9m = 30°, averaged over 100 samples, and their 
integrals normalized to unity. 



two well defined peaks to a single broad one. It means 
that the hyperbolic behaviour is progressively erased by 
multiple scattering. Besides, the width of the single peak 
scales like h - see figure [|. This scaling is in very good 
agreement with experiments jll||l9| , pof and rules out an 
analysis both in terms of the so-called telegrapher's equa- 
tion, which was proposed in [ p6| to interpolate between 
a locally hyperbolic behaviour and a large scale diffusive 
behaviour implied by scalar descriptions [fl5|]2^1 . 

More precisely, we show that the shape of this peak is 
not gaussian: see the inset of figure^. Rather, it is surpris- 
ingly close to the pure elastic response of a semi-infinite 
medium, which reads R|: 

2 1 

G e (x,h) = — 



irh [l + {x/h) 2 ] 2 ' 

independently of both the Young modulus and the Poisson 
ratio. Note that, because some force chains are deflected 
upwards and reach the top of the system, the integral of 
the vertical pressure over the bottom wall z = h gives a 
result slightly larger than 1 (typically 1.07 for h = 400 
and 9m = 30°). For a proper comparison with the elastic 
response, we have renormalized the data so that its in- 
tegral is unity. However, one should emphasize that the 
'backward' propagation of forces is a realistic effect. The 
top layer only remains stable as long as these forces do 
not exceed the weight of a single grain. For larger applied 
forces, one expects the top layer to become unstable. Let 
us now turn to an analytical derivation of the numerical 
results obtained above. 



3 A Boltzmann equation for force chain 
splitting 

The role of disorder on hyperbolic equations was previ- 
ously investigated in p7| , where the constant cq in the 
equation relating a zz and a xx was allowed to be randomly 
fluctuating in space. The study was restricted to small 
fluctuations (i.e. cn(r) = So + Sc(r), with 6c <C crj), where 
perturbation theory is valid. In this case, the two peak 
structure of the response function is preserved on large 
length scales, although the peaks are broadened. An un- 
controled extrapolation to strong disorder however sug- 
gested that the large scale equations might become ellip- 
tic. 

The chain splitting model proposed above allows us to 
investigate the strongly disordered situation. Our strat- 
egy is to write a Boltzmann equation for the probability 
density P(f, n, r) of finding an oriented force chain of in- 
tensity / in the direction n around the point r. A very 
important point here is that force chains can be oriented 
in reference to the boundary conditions. It is indeed a spe- 
cific property of hyperbolic equations to require boundary 
conditions either in the 'past' or in the 'future', but not 
in both (see the discussion in |]rj,|l8|). The possibility to 
orient the chains has been implicitly used in the above 
numerical scheme. 

For simplicity, we shall in the following neglect the 
chain 'merging' Y process, and furthermore that the split- 
ting is symmetric, i.e that n • ni = n • ri2 = cos 6*, so that 
fi = $2 = //2cos#. As will be clear from the treatment 
below, more complicated cases, with more than two off- 
springs or asymmetric scattering do not change the struc- 
ture of the large scale equations, provided the scattering 
is isotropic on average. (Anisotropic scattering can arise 
in granular materials with a non trivial contact texture, 
and leads to very interesting generalisations that will be 
mentioned in the conclusion). 

Assuming a uniform density of defects, the probability 
distribution P(f, n, r) obeys the following general equa- 
tion: 



dv 

P(f,n,r + ndr) = (1 - _)p(/, n , r ) 4 
2^ Jdn'Jdf'P(f,n',v)^(n',n)s(j 



r 



2cos6> 



(2) 



where A is equal to the 'mean free path' of force chains, and 
is of order l/(/j(ja ) in dimension D. The above equa- 
tion means the following: since a chain of grains can only 
transport a force parallel to itself , the direction of the 
force n also gives the local direction of the chain. Between 
r and r + n dr, the chain can either carry on undisturbed, 
or be scattered. The second term on the right hand side 
therefore gives the probability that a force chain initially 
in direction n' is scattered in direction n. This occurs with 
a probability ^-<^(n',n), where tf' is the scattering cross 
section, which we will assume to depend only on the scat- 
tering angle 9, for example a uniform distribution between 
— 9m and +9m- The r)-function ensures force conservation 
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and the factor two comes from the counting of the two 
possible outgoing force chains. 

Let us now multiply Eq. (^|) by / and integrate over 
/. This leads to an equation for the local average force 
per unit volume in the direction n, that we will denote 
F(n, r). This equation reads: 



An ■V r .F(n,r 
,F(n',r 



-F(n,r) 



dn' 



n • n 



/<f(n',n) + -n-F (r), 



(3) 



where we have added the contribution of an external body 
force density F (r), and a is the size of a defect or of a 
grain. This equation is the so-called Schwarschild-Milnc 
equation for radiative transfer, decribing the evolution of 
light intensity in a turbid medium p8f . We now introduce 
some angular averages of F(n, r) that have an immediate 
physical interpretation: 



p(r) = a JdQF{n,r) 
J a (r) = a J dQ n a F(n, r) 
cr a p(r) = aD J dQ n a np F(n, r), 



(4) 



(5) 
(0) 



where Jd(2 is the normalized integral over the unit sphere, 
that is introduced for a correct interpretation of a (see Eq. 
( |TT| ) below). As will be clear from the following, J is the 
local average force chain intensity per unit surface, a is 
the stress tensor. Since n 2 = 1, one finds that Trtr = Dp, 
and therefore p is the isostatic pressure. Note that J is not 
the average local force, which is always zero in equilibrium. 
The fact that J ^ comes from the possibility of orienting 
the force chains. 



4 The hydrodynamic limit: back to 'elasticity' 

Let us indeed integrate over the unit sphere Eq. (0) after 
multiplying it by different powers of n a . Using the fact 
that tj^n^n) only depends on n • n', a direct integration 
leads to: 

AV-J = (oo-l)p, (7) 

where an is called the 'albedo' in the context of light scat- 
tering p8[, and reads: 



/ , #(n',n) , 
a = dn y '/ > !• 



(8) 



A second set of equations can be obtained by multiplying 
by n a and integrating. Using the fact that 



ai = Jdn&(ri, n) = 1, 
it is easy to show that 



g(n',n) , 
an n — = n 



n • n' 



(9) 



(10) 



Therefore, the resulting equation is nothing but the usual 
mechanical equilibrium relation: 



(11) 



This relation reflects the local balance of forces chains. 
Now we multiply Eq. (^|) by n Q n^ and again integrate. A 
priori, this introduces a new three index tensor. In order 
to close the equation, we now make the an assumption 
that is usually made in the context of light diffusion, that 
on large scales the force intensity becomes nearly isotropic 
[ p8| . In this case, it is justified to expand F(n, r) in angular 
harmonics and to keep only the first terms ^| 



aF(n,r) = p(r) + Dn ■ J(r) + .. 



(12) 



The calculation requires the following identity (note that 
the repeated indices are not summed here): 

r a pys = JdQnanpn^ns = K\ [o Q/3 o 7 ,5(l - <5 Q7 ) + 

(13) 

where K x = 1/D(2 + D) and K 2 = ZK X . The resulting 
equation finally leads: 



^Dr a j3ryg'V-yJs 



ap - a-2 
D — 1 



( Da 2 - a 



~D~' 



where 



V D ~ 1 
a 2 = (cos 2 6>)= Jdn&(ri, n) n • n' 



(14) 



(15) 



Taking the trace of Eq. (14), we find, using the fact that 



the trace of a is equal to Dp: 



XD[Ki(D - 1) + K 2 ] V - J = 

ao - «2 „ , Da 2 - do . 

F> H ^ : l)p, 



D - 1 



D - 1 



(16) 



which, using the values of K\ and K 2l can be seen to be 
identical to the first equation, Eq. (0). The final equation 
appears therefore as a 'constitutive' relation between a a [3 
and the vector J. We find: 



\D 2 K! 



and 



( Da.2—ag _ 1 | 



XD 

Dai-ag _ i 
D-l 1 - 



(V a Jf3 + V J a ), (17) 



[^(X 2 -^i)V a J Q 



DKi 



a,Q — a 2 



(D-l)(ao-l) 



V- J 



(18) 



2 Actually, the third term in the expansion, proportional to 
the stress tensor, does not contribute to the following closure 
equation. Only the fourth order term (in nnn) would change 
the following results. 
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Fig. 6. Poisson ratio v for a uniform distribution of scattering 
angles, as a function of the maximum angle 9m, both for D — 2 
(plain line) and D = 3 (dashed line). 



(here a is not summed over.) These equations have ex- 
actly the canonical form of elasticity theory, provided one 
identifies the vector J with the local displacement, up to a 
multiplicative factor. This free multiplicative factor means 
that the effective Young modulus cannot be uniquely de- 
fined, as expected since it is not a dimensionless quantity. 
However, the Poisson ratio v can be defined, and reads: 

1 D - 1 + 3a - ( D + 2)a 2 
D-lD + l + a -{D + 2)a 2 ' 1 ' 

Note that this quantity is purely geometric, and does not 
involve the scattering mean free path A. Since ao > 1 and 
0-2 < 1, it is easy to show that v > 1/(D — 1) in all di- 
mensions, therefore violating the thermodynamical bound 
v <1/{D — 1) that holds for conventional elastic bodies. 
This means that the stress equations cannot be derived 
from the minimisation of a free energy functional. There- 
fore, granular materials (if describable by the present the- 
ory) should have an anomalously high Poisson ratio. The 
Poisson ratio for a uniform distribution of scattering an- 
gles is plotted as a function of the maximum angle 9 m in 
figure H, both for D — 2 and D = 3. In the weak scattering 
limit (& M -> 0), one finds v = (D + 5)/{D + 3)(D - I) 
independently of the precise shape of 'P. In the strong 
scattering limit where 9m — > tt/2, outgoing force chains 
can carry arbitrarily large forces and therefore ao — > oo. 
In this case, v —* 3/(1? — I). Such anomalous values of 
the Poisson ratio in granular materials were discussed in 
recent papers |29| , although the underlying mechanisms 
are different from the one proposed here. 

Since the above equations are formally identical to 
those of classical elasticity, one can show that V 2 p = 0, 
and V 4 J = 00. One can therefore compute the response 
function G(r) to a localized force at r = in the z di- 
rection, which is given in D = 2 by Eq. (]!]), and is, as 



figure shows, in extremely good agreement with the nu- 
merical simulation of our force chain splitting model. This 
response function is also very close to the one measured ex- 
perimentally ^,|l9|,[l6) or numerically in disordered pack- 
ingsf2l. 



5 Conclusion and extensions 

In this paper, we have shown that in the presence of strong 
disorder that 'scatters' the force chains in an assembly 
of undeformable, cohesionlcss grains, the local hyperbolic 
character of the stress equations is unstable at large length 
scales, where the equations are akin to those of a conven- 
tional elastic body. This is at variance with the weak dis- 
order calculations presented in p7| , and shows that the 
unintuitive two peak nature of the response function is 
perhaps restricted to very ordered packings, or friction- 
less particles, for which A — > oo. 

Let us stress once more that our result is not trivial 
because no displacement field is introduced in the above 
derivation (nor in our numerical model). The basic as- 
sumption of our model is the existence of local force chains, 
which have a well defined identity over several grain sizes 
a, such that this is the condition under which the 

above Boltzmann description of the force chain scattering 
is justified. Interestingly, these chains that propagate the 
forces parallel to themselves are the signature of a locally 
hyperbolic equation, and are associated with its charac- 
teristics [jl2],p . 

This work can be extended in several directions. First, 
and most importantly, the extension to non-isotropic scat- 
tering should be worked out. This means that the scat- 
tering function ^(n, n') could depend not only on n • n' 
but also on N • n', where N is a (space dependent) vec- 
tor that describes the local texture, i.e. the non-isotropic 
distribution of contacts between grains. In this case, com- 
binations like N a Jj3 or N a Np can appear in the stress 
tensor a a p. The vector N is expected to encode the con- 
struction history of the material, and should be non triv- 
ial, for example, in a sandpile constructed from a point 
source. The modification of the above theory to account 
for texture effects is needed to understand the existence of 
the pressure dip underneath the apex of the pile lIT^JTif , 
which cannot be reproduced by an isotropic elastic theory. 
In the presence of a long-ranged texture field N, the re- 
sponse function will depart from that of elasticity theory. 
Second, the fluctuations of stresses in granular materials 
are very important and have been discussed theoretically 
and experimentally [fl5| . The above theoretical framework 
should allow to make some progress, by studying higher 
moments of the local force distribution P(f, n, r), beyond 
the average force F(n, r) studied here. Third, the pseudo- 
elastic theory established above is only expected to hold 
on length scales such that L 3> A. This does not necessar- 
ily hold for small piles or small silos, where the hyperbolic 
approach seems to give very good results 0. It would 
be interesting to work out in more details the crossover 
regime where L/X is not large compared to one, using for 
example the formalism developed in |pfj. Finally, let us 
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note that the above theory should describe the large scale 
properties of random networks of rigidly jointed struts. 
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